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ABSTRACT 

The  problem  of  linearly  constrained  least  squares  has  many  applications  in  signal  processing.  In  this  paper, 
we  present  a  perturbation  analysis  of  a  linearly  constrained  least  squares  algorithm  for  adaptive  beamforming.  The 
I>erturbatiua  bounds  for  the  solution  as  well  as  for  the  latest  residual  element  are  derived.  We  also  propose  an  error 
estimation  scheme  for  the  residual  element,  which  can  be  incorporated  into  a  systolic  array  implementation  of  the 
algorithm. 


1.  INTRODUCTION 

The  least  squares  problem  with  linear  equality  constraints  has  important  applications  in  signal  processing, 
e.g.,  adaptive  beamforming.  To  solve  this  problem,  MeWhirter  and  Shepherd  [5]  proposed  a  systolic  algorithm  and 
architecture.  In  this  paper,  we  present  a  perturbation  analysis  of  the  problem  and  propose  an  error  estimation 
scheme  for  the  MeWhirter-Shepherd  (MS)  algorithm  [5).  This  paper  is  organized  as  follows.  The  least  squares 
problem  is  defined  in  Section  2  and  error  bounds  are  derived  in  Section  3.  An  error  estimation  algorithm  is  given 
in  Section  4,  and  in  Section  5  a  numerical  example  is  presented  to  illustrate  how  well  our  new  algorithm  works. 

2.  PROBLEM  DEFINITION 

Given  an  n  x  q  complex  data  matrix  A'(n),  the  least  squares  problem  with  linear  equality  constraints  is  to  find 
a  q-element  complex  vector  w(n)  such  that 

llA(n)u;(n)l|  =  min  (2.1a) 


subject  to  the  linear  constraints 


Sw{n)  =  6, 


(2.1b) 


where  S  is  &  k  x  q  {k  <  q)  complex  matrix  and  6  is  a  Ir-element  complex  vector.  Throughout  this  paper,  we  use 
the  2-norm: 


II  11  =  11  lb- 


In  signal  processing,  new  data  arrives  continuously.  Define  the  data  matrix  X{n)  recursively  by 


i.e.,  the  nth  row  x(n)^  represents  a  snapshot  at  time  n.  Our  goal  is  to  compute  the  n-th  residual  element 

r„  =  x(n)^u<(n).  (2.2) 
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Is  the  solution  vector  w(n)  unique?  Define  a  (fc  +  n)  x  5  matrix  Sx{n)  by 


Sx{n)  -  (  Y(n))  • 

We  assume  that  k  +  n  >  q.  The  solution  is  unique  if  and  only  if  the  matrix  Sx{n)  has  full  column  rank;  that  is, 
the  overdetermined  matrix  equation 

Sx(n)w{n)  =  0  (2.3) 

has  a  unique  solution  w(n)  =  0. 

Next,  we  wish  to  transform  (2.1)  into  a  familiar  unconstrained  problem;  see  [3]  and  [4].  Let 

p  =  q-k 


and  partition  the  matrix  S  ils 


S  =  (Si  S2). 


where  Si  is  k  x  k  and  Sn  is  k  x  p.  For  simplicity,  we  assume  that  Si  is  nonsingular  and  upper  triangular;  for 
example.  Si  may  be  the  result  of  an  initial  QR  decomposition  of  S.  Accordingly,  we  also  partition  X(n)  as 

A'(n)  =  ( ,Yi(/i)  A'2(n)), 

so  that  A’l  is  n  X  k  and  X2  is  n  x  p.  Then  (2.3)  becomes 


which  is  equivalent  to 


where 


Af(n)  Af(n))*"^"^-°’ 
) 


5i  S2 
0  C(n) 


ti;(n)  =  0, 
C(n)  =  A'2(n)  —  A'i(n)Sf ‘S2. 


The  matrix  C(tj)  is  called  the  Schur  complement  of  Si  in  Sx-  The  equation  (2.3)  has  the  trivial  solution  if  and 
only  if  C(n)  has  full  column  rank.  We  proceed  to  eliminate  the  constraints.  Let 


so  that  tei(n)  is  ib  x  1  and  u;2(n)  is  p  x  1.  Since 

Si«/i(n)  +  S2«'2(«)  = 


we  get 

u;i(n)  =  Sf ‘6  -  Sf  ^S2U^2(«)- 

(2.4) 

Let 

v(n)  =  -Ai(n)Sf ‘6. 

We  derive 

llC(n)u;2(n)  -  v{n)\\  =  min, 

(2.5) 

an  unconstrained  problem  analyzed 

in  [3],  [4].  Now,  what  about  the  residual  element  r„? 

Define  the  Schur 

complement  matrix  C(n)  recursively  by 


_  /C(n-1)\ 


Partition  the  row  vector  x(n)^  so  that 


x{nf  =(xi{nf  Z2{n)^), 
where  Xi(n)^  is  1  x  ife  and  x^{n)^  is  1  x  p.  We  get 

c(n)^  =  X2{nf  -  Xi(n)^Sf ‘52- 

Let  v„  denote  the  n-th  element  of  v(n).  The  last  residual  element  of  (2.5)  is  then 

c(n)^u;2(n)  -  d„  =  x^in)'^ W2{n)  +  xiirifwiin)  +  v„-v„  =  r„, 

i.e.,  the  same  residual  element  as  desired  by  the  constrained  problem  (2.1). 

How  do  we  calculate  r„  recursively?  Suppose  that  we  have  available  a  QR  decomposition  of  the  (n  —  1)  x  p 
matrix  C(n  —  1): 

C{n  -  1)  =  Q(n  —  l)R{n  —  1), 

where  Q(n  —  1)  is  (n  —  1)  x  p  with  orthonormal  columns  and  the  matrix  R{n  —  1)  is  p  x  p  upper  triangular.  The 
problem  (2.5)  is  reduced  to 

where  u(n  —  1)  =  Q(n  —  l)^v{n  —  1).  We  triangularize  the  coefficient  matrix  by  a  unitary  matrix  P.  Then 

HfR{n-l)  u{n  -  l)\  _  f  R{n)  u(n)\ 

I,  c(n)^  v„  )  0^  7  /  ’ 

so  that  R(n)  is  p  x  p  upper  triangular.  The  matrix  P  consists  of  p  Givens  matrices.  From  P  and  Q(n  —  1)  we  can 
construct  an  n  x  p  orthonormal  matrix  Q(n)  such  that  C(n)  =  Qln)R{n)  and  u(n)  =  Q(n)^ v{n).  The  desired 
element  Vn  is  given  by 

rn  =  -(ci  ...Cp)7, 

where  ci , . . . ,  Cp  denote  cosines  of  the  p  rotations  that  make  up  P. 

3.  PERTURBATION  ANALYSIS 

Eld4n  [1]  presented  a  perturbation  analysis  of  the  linearly  constrained  leaist  squares  problem.  Since  his  theory 
is  general,  it  involves  weighted  pseudoinverses  and  their  corresponding  condition  numbers.  In  this  section,  we 
derive  simpler  perturbation  bounds  for  the  solution  w{n)  as  well  as  for  the  residual  element  r„.  To  simplify  our 
presentation,  we  will  drop  the  argument  (n)  for  the  matrices  and  vectors,  and  let  /c(M)  denote  the  condition  number 
of  a  matrix  M  with  respect  to  the  2-norm. 

Let  w  solve  the  perturbed  least  squares  problem 

II  (  Xi -f- f  L’x,  X2  +  )  ihll  =  min  (3.1a) 

subject  to  the  perturbed  linear  equality  constraints 

(Si+(Esi  S2  +  (Es,)w  =  b  +  tfi,.  (3.1b) 

Suppose  that  t  >  0  is  a  real  variable  and  let 

C  +  tEc  =  (X2  -(-  tEx,)  -  (Xi  +  tEx,)(Si  -I-  tEs,)-HS2  +  tEs,) 

and 

V  +  tft,  =  — (Xi  +  tEx,  ){Si  +  tEs^  +  tft). 


3 


Recall  that  Si  is  nonsingular  and  that  C  has  full  column  rank.  Suppose  e  is  sufficiently  small  so  that  for  t  €  [0,f] 
we  have  Si  +  tEg,  is  nonsingular  and  C  +  tEc  has  full  column  rank.  Let  w(t)  solve  the  matrix  equation 


/Si  +  tEs,  Si+tEs^  ^  +  A 

0  +  +  \(C  +  tEcf(v  +  tfj)- 


(3.2) 


Then  u;(0)  and  w{€)  are  solutions  to  problems  (2.1)  and  (3.1),  respectively.  Define  w  =  u;(0)  and  it’  =  U'(e).  Then 

w  =  u>(0)  +  fu)(0)  +  0(e~). 

Differentiate  (3.2)  with  respect  to  t  and  set  t  =  0.  We  get 


Let 


Then 


(  o‘  +  C"c)  {e^v  +  C^/, 


(‘o'  C%)  '=S-\C"Cr^  and  ||C||<||C||. 

Solving  for  w{Q)  in  (3.3),  we  obtain 

/  c  ^  \  ^  r 

ti^(o)  =  ( ; 

=  S-'(C'"C)-'(?" 

where  r  =  Cw^  —  v  denotes  the  residual  vector.  Furthermore,  by  assuming 

IIM1<IH1.  I|£:cll  <  licil 

and 

i(t  toil -1(^0  ?)||<II^INIC||, 

we  derive  the  inequality 

N(0)||  <  115-^11  ||(C'^C)-‘C"||  [||d||  +  !1(?|1  IISIl  ||u,|l]  +||S-'||  ||(C"C)-‘||  ilCII  ||r||. 

Consequently,  we  obtain  the  following  perturbation  result. 

Lemma.  Using  the  notations  defined  above  and  assuming  that  (  in  (3.1)  is  sufficiently  small  so  that  the  inequa/ities 
(3.5)  are  satisfied,  we  get 


0  0 

0  Efic  ' 


(33) 


(3.4) 


(3.5a) 


(3.5b) 


llih- 


tn 


<  c  |/c(5)«(C)  (  .  ^^3 - +  l)  +  «(S)«(C)^-^— - \ 

-I  viiqilisiiiHi  J  liqiPiiikiii 


+  0(€^).a 


(3.6) 


To  illustrate  the  effect  of  k:(S)  on  the  solution  of  (2.1),  consider  a  simple  example  in  which  S  =  (  Si  0  )  and 
X  =  ( I„  In),  where  /„  is  an  n  x  n  identity  matrix  and  n  <k  <2n.  By  observation,  mi  =  Sf*6  and  mo  =  -U’l. 
Since  k(S)  =  «(Si)  in  this  example,  we  see  why  the  presence  of  ic(S)  is  necessary  in  (3.6). 
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We  proceed  to  derive  a  bound  for  the  error  in  the  residual.  Let 


\mj  V  0 


S2 + ^£53 

C+tEc  ) 


U)(0 


differentiate  the  equation,  and  then  set  <  =  0.  Using  (3.4)  to  substitute  for  ii;(0),  we  get 

(4))=("»' 


=  (I-CC^) 


Esi  Esj 

0  Ec 


y- 


u 


H  /^\- 1 


-C(C^C) 


0 

E}ir 


where  =  (C^ C)  .  Consequently, 

r(0)  =  (/  -  CCy(EcW2  -  -  CiC^Cy^E^r. 

As  for  the  residual  element  we  have  r„  =  e^r,  where  e„  =  {0,...,0,1)  denotes  the  n-th  unit  coordinate  vector. 
Using  the  assumptions  (3.5a)  smd  noticing  that  W2  =  C^v  and  »•  =  (/  —  CC^)j;,  we  derive  our  major  result. 

(3.7) 

(3.8) 


Theorem.  Under  the  same  conditions  as  in  the  Lemma,  we  get 

Ilf'  -  r|| 


IIHI 


<€[(|/-CCMl(2K(C)  +  l)]+0(f2) 


and 


|f‘n  -  '•nl 


<  €  [||/  -  CC»||(K(C)  +  \\C\\  l|C»e„||  +  1)]  +  0{e^).  a 


Here  are  some  additional  remarks.  If  we  set  5  =  /  and  A  =  0,  then  (3.6)  leads  to  a  perturbation  bound  for  the 
standard  least  squares  problem  [2],  We  also  note  that  j|C5u;||*  +  |lr||^  =  |ld|p.  Thus,  we  can  define 

cos0  =  i|C5u>||/lld|| 

and  use  (1/cosfl)  and  tan0  in  (3.6).  The  bound  (3.7)  is  similar  to  a  result  derived  in  [2].  The  inequality  (3.8) 
indicates  that  l^n  —  r„|  depends  on  /c(C)  as  well  as  on  |ln||.  Both  (3.7)  and  (3.8)  can  be  simplified  by  using  the 
relation  that  ((/  —  CC^||  =  min{l,  n  —  p}. 


4.  ERROR  ESTIMATION 

Although  the  error  bound  (3.8)  is  simple,  it  requires  C^e„,  whose  computation  involves  at  least  a  back-solve. 
In  this  section,  we  present  an  error  estimation  scheme  for  the  desired  residual  element.  When  the  new  data  vector 
x(n)^  arrives,  it  is  first  processed  by  5  so  that  Xi(n)^  is  annihilated.  In  particular,  let 


(0) 

1 


rj°*)  =  x(n)^  and  u*°’  =  0. 


Then  the  preprocessing  proceeds  as  follows: 


Sl.l  Sii  +  i 

.  ‘1+1 


s/.f')  (  \  q\  f  si,i 


S|.(+l 
)  ,(<-!) 
‘1+1 


,U-1) 


and 


b, 

u(') 


)^(4,  nU-)' 
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for  I  =  1,2, . . . ,  fc,  where  gi  =  zj*  Writing  in  algorithmic  form,  we  have 

for  /  =  1, 2,  . .  ,k 


9i  -  : 

for  j  =  /  +  1, . . .  ,9 


-  9,SU 


=  «('  —  gibi 


The  above  process  shows  that 

(4+1  •••  •g'"’)  =  ■^2(n)^  -  xi(n)^S;''S2  and  u‘*’ = -Ji(n)^5f  ^6. 

These  two  variables  are  then  used  for  updating  the  QJi  decomposition  of  C(n  —  1)  and  computing  the  residual 
element.  We  present  below  the  algorithm  derived  in  [4]. 

for  /  =  1, 2, . . .  ,p 


cos^i  =  ! 

sin^/  =  ; 

for  j  =  /  +  1 , . . . ,  p 


c5"’  =  c{"  cos  6i  +  *'  sin^j ; 

(k+l)  (n-l)  •  a  1  (t+l-1)  a 

=  -c)j  sm  Bi  +  2^+^  cos 


(t+i-i). 


=  nj"  cos^i  +  ^*sin0|; 

j,(*+0  _  sin^i  +  «^‘''''‘~^^cos0( 


r„  =  u<*+P)n^^,cos0,. 

In  the  above,  cj*'  (for  ifc  =  n  —  l,n)  denotes  the  (i,  j)-element  of  C{k)  and  the  j-th  element  of  v(k). 
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Now,  we  discuss  an  error  estimation  scheme  for  the  preprocessing.  Let  '  denote  the  corresponding  computed 
value  and  //  the  floating  point  computation.  In  the  above  procedure  we  calculate 

91  = 

Define  the  relations  between  the  exact  and  computed  quantities  as  follows: 

91  =  <;/(!  +aj6(f)). 

6/  =  6/(1  +/i/6/(e)), 

where  |<At,j(f)|  =  0(e),  =  0(e),  |4/(f)|  =  0(e),  =  0(e)  and  |6;(f)|  =  0(e).  The  five  quantities  (t,  j, 

Cj*\  a/,  and  fii  are  all  real  and  nonnegative.  We  also  assume  that  the  errors  such  as  <Ti,j<t>i,j(e)  and  (e) 

are  so  small  that  higher  order  terms  like  (<rtj<t>ij(e))-  and  (ffij<t>ij(e))((ll‘~^\>\‘~^\e))  can  be  ignored.  Using  the 
lemma  in  [3],  we  obtain  the  following  algorithm  for  estimating  the  errors  in  preprocessing. 

for  i  =  1,2, ...  ,k 


ai  =  max{C/  "^'.tr/,/}  ; 
for  j  =  /  +  1, . . . 


As  explained  in  [3],  the  above  estimation  scheme  can  be  incorporated  with  the  preprocessing  procedure  and  imple¬ 
mented  on  the  same  systolic  architecture.  Additional  time  is  minimal  because  the  calculations  can  be  carried  out 
during  the  otherwise  idle  time  of  the  processors. 

The  error  estimate  for  r„  can  be  obtained  by  the  algorithm  presented  in  [3]  using  (C^+i  •  ■  *)  and  q'*’  as 

the  error  estimates  for  •••  ■2'^^)  and  respectively.  Again,  we  list  the  error  estimation  algorithm  and 

refer  the  details  to  [3].  Define  the  relations  between  the  exact  and  computed  quantities  as  follows: 


+Ci,p+ia,,p+i(f)), 


cosdi  =  COS0/(1  +  <T/,/(6[y(f )), 
siuBi  =  sin0/(l -I- cr/,;0}y(e)), 

v(  ^  *(1  -f  (T,  p+i(^/  )), 


f‘n  =  2-„(l  q^(f)). 
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The  following  algorithm  estimates  the  error  in  the  last  element  of  the  residual  vector: 


for  /  =  1, 2, . . .  ,p 
begin 


erij  =  max{^,,,,(;itV~^’}; 

for  j  =  /  +  1 , . .  . ,  p 


begin 


|cj“~  '*  cot  Si  max{(i 
.1  =  — ^ - TT 


n+i2<*+'- 


, •(!  +  *) 
'’t+j 


’  sin«i  tnax{0.j,ff|.i)|+| 


.(K+'-i)  . 


end; 


_  maxf^i  ,.).|,g|.i}|+|u‘*~'’'~‘*sin»i  maxlt)***'" ‘\g|,i )| 

~  co« sin 

(1+/)  _  Iti}**"**  sin»[  max((i,,.(.|,<7i.i)|4-|u**'''*~**  cog»i  max{i)**~*’‘~ '* ,g|,i}| 
~  ~  ain  ®i  — «**+*“  **  co«  #r| 

end; 

T)  =  »7(*+P)max{<Ti,i,  •  •  .,(Tp,p}  . 


5.  AN  EXAMPLE 

The  example  in  this  section  shows  that  the  computed  residual  element  may  be  accurate  even  when  the  matrix 
C  is  Unconditioned.  In  this  case,  the  proposed  scheme  gives  a  better  error  estimate  than  (3.8).  Both  the  MS 
algorithm  and  the  error  estimation  algorithm  were  implemented  using  MATLAB  and  run  on  a  VAX  8300  with 
machine  precision  f  =  1.1102  x  10"'®  in  the  Communications  Research  Laboratory  at  McMaster  University. 

Example.  Suppose  the  exact  constr^nt  matrix  and  corresponding  right  side  vector  are 

(10^  0  0  1  0  0\  / -120000  v/2/7\ 

0  10"^  0  0  1  0  j  and  6  =  I  v'iO/700 

0  0  1  0  0  1/  \  6x/5/7  / 

Thus  we  set  the  error  estimates  as  (Ti  j  =  pi  =  1  and  d, ■,;(€)  =  6i{()  =  e.  The  data  matrix  at  time  n  —  1  is 

(-1  -x/5  -2\/r0  0  0  0\ 

0  -1  n/2  0  0  0  j  . 

0  0  -1  0  0  0/ 

Suppose  we  know  the  exact  R(n  —  1)  and  u(n  —  1): 


/  0.001 

1000\/5 

2\/l0  \ 

/-10v/2/7\ 

ft(n-  1)  =  0 

1000 

-2V2 

and  u{n  —  1)  =  j  4\/l0/7  j 

\  0 

0 

1  / 

\  6v^/7  / 

Similarly,  the  error  estimates  of  their  elements  are  all  initialized  as  e.  Now  the  new  data 

x(n)^=(-l  -Vb  -2v/l0  0.001  0  0  )  and  =  0 
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are  available  and  their  error  estimates  are  initialized  as  =  1,  for  j  =  1 . 6,  and  r;'®*  =  1,  respectively.  After 

preprocessing,  we  get  c(n)^  =  (0.002  lOOOv/S  2\/To)  and  v„  =  — lOv^/7.  The  corresponding  error  estimates 
are  =  1,  for  j  —  4,5,6,  and  =  23.  The  QR  updating  scheme  and  its  error  estimation  algorithm  are  then 

applied  to  R{n  —  1),  u(n  —  1),  c(n),  ii„  and  their  error  estimates.  The  exact  residual  element  r„  =  6\/2/35.  The 
computed  error  is 

!»>„  -r„l  =  1.11  X  10-'®. 

The  condition  number  of  C(n)  is  4.6  x  10®  and  the  error  bound  as  given  by  (3.8)  equals  3.40  x  lO”'*.  The  estimation 
algorithm  gives  a  much  more  accurate  value  of  9.62  x  10"'®. 

ACKNOWLEDGEMENTS 

The  work  of  F.  T.  Luk  was  supported  in  part  by  the  Rensselaer  Polytechnic  Institute,  and  by  the  Army  Research 
Office  under  grant  OAAL03>90-G-0104  and  the  Joint  Services  Electronics  Program  under  contract  F49620-90-C- 
0039  at  Cornell  University.  The  work  of  S.  Qiao  was  supported  in  part  by  Natural  Sciences  and  Engineering 
Research  Council  of  Canada  under  grant  OGP0046301. 

5.  REFERENCES 

[1]  L.  Eld^n,  “Perturbation  theory  for  the  least  squares  problem  with  linear  equality  constraints,"  Si.\M  J.  Numer. 

Anal.,  Vol.  17,  No.  3,  June  1980,  pp.  338-350. 

[2]  G.  H.  Golub  and  C.  F.  Van  Loan,  Matrix  Computations,  Second  Edition,  The  Johns  Hopkins  University  Press, 

Baltimore,  Maryland,  1989. 

[3]  F.  T.  Luk  and  S.  Qiao,  “Analysis  of  a  recursive  least-squares  signal-processing  algorithm,”  SIAM  J.  Sci.  Stat. 

Comput.  Vol.  10,  No.  3.  May  1989,  pp.  407-418. 

[4]  J.  G.  McWhirter,  “Recursive  least-squares  minimization  using  a  systolic  array,”  Real  Time  Signal  Processing 

VI,  K.  Bromley,  Ed.,  Proc.  SPIE  431,  1983.  pp.  105-112. 

[5]  J.  G.  McWhirter  and  T.  J.  Shepherd.  “A  systolic  array  for  linearly  constrained  least-squares  problems.” 

Advanced  Algorithms  and  Architectures  for  Signal  Processing,  J.  M.  Speiser.  Ed.,  Proc.  SPIE  696.  1987.  pp.  80- 

87. 


9 


